#ifndef MATRIX_H
#define	MATRIX_H

#ifdef	__cplusplus
extern "C" {
#endif
#define matrix(m,i,j,n) m[(n)*(i)+(j)]

void matrixMul(int* a, int* b, int* result, int size) { //basic matrix multiplication - FFT would indeed be better
        int* buf = (int*) malloc(sizeof (int) *size * size);
        int i, j, k;
        for (i = 0; i < size; i++) {
            for (j = 0; j < size; j++) {
                matrix(buf, i, j, size) = 0;
                for (k = 0; k < size; k++) {
                    matrix(buf, i, j, size) += matrix(a, i, k, size) * matrix(b, k, j, size);
                }
            }
        }
        for (i = 0; i < size * size; i++) {
            result[i] = buf[i];
        }


        free(buf);
        //printf("done matrix mul \n");
    }

    void matrixPow(int *a, int k, int*result, int size) { // determination of the number of directed walks from one node to the other by adjacency matrix potentiation
        //printf("start matrix pow \n");
        int j = 0;
        int i, exp;
        for (i = 0; i < size * size; i++) {
            result[i] = a[i];
        }
        exp = k - 1;
        while (exp != 0) {
            if ((exp % 2) == 1) {
                matrixMul(a, result, result, size);
                exp--;
            }
            matrixMul(a, a, a, size);
            exp = exp / 2;
        }
        //printf("done matrix pow \n");
    }
    
#ifdef	__cplusplus
}
#endif

#endif	/* TOOLS_H */